Tomographic images are produced by converting observed projections (data) into an image. For example, in x-ray CT imaging, x-ray beams are directed at an object and the beams are attenuated by various amounts due to varying structures within the object. On the other side of the object, the attenuated beams are measured with detectors. Such projections are produced at many different angles around the object. Not only are these measurements noisy, the relative noise level depends upon the amount of attenuation. Projections through dense materials such as bone and especially metal have lower signal-to-noise ratios than projections through flesh, water, or other less dense materials. Coping with the large and spatially varying fluctuations in the number of detected photons often requires a statistical smoothing technique, also known as regularization, to improve the image.
In the absence of noise, a basic problem of CT scans amounts to determining an unknown image f=f(x, y) from its projections, or Radon transform Rf where (.Rf)(t,θ):=∫∫f(xy)δ(t−x cos θ−y sin θ) dxdy, where (x, y) are planar coordinates, t is the location along each projection, and θε[0, π] is the orientation of the projection. Unfortunately, because measurements are never perfect, what is actually observed is the noisy projection data g=g(t,θ). To emphasize the essentials of the tomography problem, the unknown f=f(x,y) and the observed g=g(t,θ) are viewed as functions, although the implementation is discrete. Next, assume the standard independence and locality assumptions that the likelihood P(g|f) or conditional distribution of g given f, equals Πp(g(t,θ) |(Rf)(t θ,)), where the product is over all (t,θ). This product may be nonzero if a finite number of factors is taken. Here, the observed number of x-ray photons b exp(−g(t,θ)) is Poisson distributed with mean w(t,θ):=b exp(−(Rf)(t,θ)), where b is the mean number of incident photons and f(x, y) is the attenuation coefficient at (x, y). The task is to infer the image f given noisy projections g.
Because this inverse problem is ill-posed, one typically imposes extra constraints on f. In regularization, a type of penalized maximum likelihood, inferring f amounts to finding that f which minimizes−ln P(g |f)+p(f), where p(f) characterizes the extra constraint on f Here we impose smoothness via the gradient
      ∇    f    =      (                            ∂          f                          ∂          x                    ,                        ∂          f                          ∂          y                      )  with the quadratic penalty p(f)=β∥∇f∥2, where β>0 and ∥f∥2=(f,f), using the inner product (f1,f2)=∫f1f2. Following Saver et al., “A local update strategy for iterative reconstruction from projections”, IEEE Transactions of signal Processing, vol. 41, no. 2, pp. 534–548 (1993), and to simplify the presentation, I approximate×ln P(g|f) with the quadratic form ∥g−Rf∥w2 where ∥g∥w2=(g, Wg) is a weighted norm with (diagonal) weight operator W satisfying (Wg)(t,θ)=w(t,θ)g(t,θ). The weight w(t,θ) is small for those rays passing through dense materials such as bone or metal, and larger otherwise. This formulation of tomography requires solving the following “hard” optimization problem.Problem 1 (Coupled Regularization). Given the projection data g, find the image f that minimizes∥g−Rf∥w2+β∥∇f∥2.(Because w depends on Rf, which is unknown, we use g instead to obtain a quadratic functional above. In general, however, one could solve the full nonlinear problem and only make this substitution as a first guess of w. One can compute f iteratively, using the current best estimate of Rf to set w.) To proceed, let Δ denote
                    ∂        2                    ∂                  x          2                      +                  ∂        2                    ∂                  y          2                      ,the Laplacian in the plane. Using integration by parts and zero boundary conditions, recall that ∥∇f∥2=(f,−Δf). Then the (linear) Euler-Lagrange equation for Problem 1 is
                                                                        R                *                            ⁢              WRf                        -                                          β                ⁡                                  (                                                                                    ∂                        2                                                                    ∂                                                  x                          2                                                                                      +                                                                  ∂                        2                                                                    ∂                                                  y                          2                                                                                                      )                                            ⁢              f                                =                                    R              *                        ⁢            Wg                          ,                            (        1        )            where f is unknown and A* denotes the adjoint (or transpose) of linear operator A (R* is also known as the backprojection operator). By examining equation (1), we see that Problem 1 is hard in two related ways.
First, the problem constraints occur in two different domains. Fidelity to the data (∥g−Rf∥w2) is enforced in the Radon domain {(t,θ)}, while smoothness (∥□f∥2) is imposed in the image domain {(x, y)}. Thus in equation (1) the operators R and R* are for shuffling back and forth between these domains; iterative solution techniques typically compute these forward and backprojections explicitly and often at great expense.
Second, observe that equation (1) is a coupled equation in the two-variable functions, i.e., in the large set of variables {f(x,y), for all x,y} under some reasonable discretization of x and y. The coupling arises first because both x and y derivatives are present; in addition, R and R* are integral operators, and so are not even local. The computational difficulty in solving equation (1) and related tomographic problems (e.g., emission) has spawned a great deal of work in optimization.
Thus, the need exists for a processing technique which can be formed in a single domain.